The infrared side of galaxy formation. 
I. The local universe in the semi- analytical framework. 
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ABSTRACT 



We present a new evolutionary model for predicting the far-UV to sub-mm properties of the 
galaxy population. This combines a semi-analytic galaxy formation model based on hierarchi- 
cal clustering (GALFORM, Cole et al. 2000) with a spectro-photometric code which includes 
dust reprocessing (GRASIL, Silva et al. 1998). The former provides the star formation and 
[ metal enrichment histories, together with the gas mass and various geometrical parameters, for 

■ a representative sample of galaxies formed in different density environments. These quantities, 

together with a few other assumptions concerning the spatial distribution of dust and its optical 
properties, allow us to model the spectral energy distributions (SEDs) of galaxies, taking into 
account stellar emission and also dust extinction (absorption plus scattering) and re-emission. 
In the spectro-photometric code dust is considered only in the disk, but the general radiation 
field is contributed by both the disk and the bulge components with their own distinct age and 
mctallicity distributions. Two phases are considered for the dust: molecular cloud complexes, 
where stars are assumed to be born, and the diffuse interstellar medium. The model includes 
both galaxies forming stars quiescently in disks, and starbursts triggered by galaxy mergers. We 
A | test our models against the observed spectro-photometric properties of galaxies in the local Uni- 

verse, assuming a CDM cosmology with Qq = 0.3 and Ao = 0.7. The models reproduce fairly well 
\ the SEDs of normal spirals and starbursts from the far-UV to the sub-mm, and their internal 

extinction properties. The starbursts follow the observed relationship between the FIR to UV 
luminosity ratio and the slope of the UV continuum. They also reproduce the observed starburst 
attenuation law (Calzetti et al. 1999). This result is remarkable, because we use a dust mixture 
which reproduces the Milky Way extinction law. It suggests that the observed attenuation law 
is closely related to the geometry of the stars and dust. We compute galaxy luminosity functions 
over a wide range of wavelengths, which turn out to be in good agreement with observational 
data in the UV (2000A), in the B and K bands, and in the IR (12 - 100/mi). Finally, we inves- 
tigate the reliability of some star formation indicators which are based on the properties of the 
continuum SEDs of galaxies. The UV continuum turns out to be a poor star formation indicator 
for our models, whilst the infrared luminosity is much more reliable. 

Subject headings: Galaxies: evolution; Galaxies: formation; Galaxies: fundamental parameters; Galaxies: 
interactions; Galaxies: starburst; Infrared: galaxies; Ultraviolet: galaxies 
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different wavelengths out to z < 5. By combining 
observations in the UV, optical, IR and sub-mm, 
we can now start to reconstruct the history of 
star formation in galaxies over the epochs when 
the bulk of the stars have formed (e.g. Madau 
ct al. 1996; Steidel et al. 1999; Hughes et al. 

1998) . On the theoretical side, models based on 
the paradigm of structure formation through hi- 
erarchical clustering (which has successfully con- 
fronted a wide range of observations on large scale 
structure and microwave background anisotropies) 
have now been developed to the point where they 
can make definite predictions for the observable 
properties of galaxies (luminosities, colours, sizes, 
morphologies etc) at all redshifts, starting from an 
assumed initial spectrum of density fluctuations. 
The key technique for making these predictions 
has been that of semi- analytical modelling (White 
& Frenk 1991; Lacey & Silk 1991; Kauffmann et 
al. 1993; Cole et al. 1994; Somerville & Primack 

1999) . In this technique, one applies simplified an- 
alytical descriptions of the main physical processes 
of gas cooling and collapse, star formation, feed- 
back effects from supernovae, galaxy merging etc, 
with the backbone being a Monte Carlo descrip- 
tion of the process of formation and merging of 
dark matter halos through hierarchical clustering. 
The predicted star formation histories are then 
combined with detailed stellar population models 
to calculate galaxy luminosities at different wave- 
lengths. Conversely, direct numerical simulations 
have been enormously successful in studying the 
evolution of structure in the dark matter on a huge 
range of scales (e.g. Jenkins et al. 1998), but cur- 
rently do not have sufficient spatial resolution to 
simultaneously follow all the processes involved in 
galaxy formation. 

The semi-analytical models have been success- 
ful in predicting and/or explaining a large range of 
galaxy properties, both at low and high redshift, 
for instance, luminosity functions and colours in 
different optical and near-IR bands (Lacey et al. 
1993; Kauffmann ct al. 1993; Cole ct al. 1994), 
the mixture of galaxy morphologies and the evolu- 
tion of elliptical galaxies (Kauffmann et al. 1993; 
Baugh et al. 1996; Kauffmann 1996), the prop- 
erties of Lyman-break galaxies at high redshift 
(Baugh et al. 1998; Governato et al. 1998), the 
sizes and circular velocities of galaxies (Cole et al. 

2000) , and galaxy clustering evolution and the na- 



ture of the clustering bias (Kauffmann ct al. 1997; 
Baugh et al. 1999; Diaferio et al. 1999; Benson et 
al. 2000). However, with very few exceptions, 
these semi-analytical models have ignored both 
extinction and emission by interstellar dust, and 
calculated only the direct stellar emission in the 
UV, optical and near-IR. This has been partly be- 
cause the importance of dust was generally under- 
appreciated, especially for high redshift galaxies, 
but also because of the lack of physically realistic 
models for predicting dust effects. 

This situation has now begun to change. On 
the one hand, there have been several observa- 
tional discoveries demonstrating the importance 
of dust effects for building a complete picture of 
galaxy formation. (1) The discovery of a cosmic 
far-IR/sub-mm background by the COBE satel- 
lite (Puget et al. 1996; Guiderdoni, et al. 1997; 
Dwek, et al. 1998; Fixsen et al. 1998; Hauser 
ct al. 1998; Schlegel et al. 1998), whose energy 
density indicates that, as suggested already by 
Wang (1991) and Franceschini et al. (1994), a 
large fraction of the energy radiated by stars over 
the history of the universe has been reprocessed 
by dust. (2) The discovery that the population of 
star forming galaxies at z ~ 2 — 4 that have been 
detected through their strong Lyman-break fea- 
tures are substantially extincted in the rest-frame 
UV (Pettini et al. 1998; Steidel et al. 1999). (3) 
The discovery of a population of sub-mm sources 
at high redshift (z > 1) using SCUBA, whose 
luminosities, if they are powered by star forma- 
tion in dust-enshrouded galaxies, imply very large 
star formation rates (~ 10 2 M Q yr _1 ), and a to- 
tal star formation density comparable to what is 
inferred from the UV luminosities of the Lyman- 
break galaxies (Smail et al. 1997; Hughes et al. 
1998; Lilly et al. 1999). (4) The ISO detection of a 
population of strong IR sources; 15 /im ISOCAM 
(e.g. Oliver, et al. 1997; Elbaz, et al. 1999) and 
175 /im ISOPHOT surveys (e.g. Kawara, et al. 
1998; Puget, et al. 1999) indicate a population 
of actively star forming galaxies at 0.4 < z < 1.3, 
which boosts the cosmic star formation density by 
a factor ~ 3 with respect to that estimated in the 
optical from the CFRS (Flores et al. 1999). For 
(1) and (3), there is the caveat that the contribu- 
tion from dust-enshrouded AGNs to the sub-mm 
counts and background is currently uncertain, but 
probably the AGNs do not dominate (Granato et 
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al. 1997; Almaini et al. 1999; Madau 1999). 
These discoveries demonstrate that in order to un- 
derstand the history of star formation in the uni- 
verse from observational data, one must have a 
unified picture that covers all wavelengths from 
the far-UV to the sub-mm. The UV and the far- 
IR are expecially important, since young stellar 
populations emit most of their radiation in the 
rest-frame UV, but a significant fraction of this is 
dust reprocessed into the rest-frame far-IR. 

On the theoretical side, it is now possible for 
the first time to construct true ab initio models in 
which the galaxy formation itself and stellar emis- 
sion and dust absorption and emission are calcu- 
lated from first principles, based on physical mod- 
els, and avoiding observational parameterizations 
for various key ingredients (e.g. shape of the lumi- 
nosity function, dependence of dust temperature 
on galaxy properties). These new models, which 
provide a unified treatment of emission from stars 
and dust, and predict the evolution of galaxy lu- 
minosities from the far-UV to the mm, are the 
subject of this paper. 

The effects of dust on galaxy luminosities at 
different wavelengths have been included in some 
previous galaxy evolution models, at various lev- 
els of sophistication, but mostly in the context 
of backwards evolution models, where one tries to 
evolve back from observed galaxy properties at 
the present day, in contrast to the semi-analytical 
models, where one evolves forward from cosmo- 
logical initial conditions. In backwards evolution 
models, one starts from the observed luminosity 
functions of different types of galaxy at the present 
day, assumes a different star formation history for 
each type, and calculates the luminosity evolution 
for each type, to predict what the galaxy popula- 
tion would have looked like in the past. Guidcr- 
doni & Rocca-Volmerange (1987) were the first 
to include dust absorption in a model of this type, 
based on a ID slab model for the star and dust 
distribution, and calculating the dust content self- 
consistently from a chemical evolution model. The 
same treatment of dust was later used in the semi- 
analytical galaxy formation models of Laccy et al. 
(1993). In both cases, the models were used to cal- 
culate galaxy luminosities and number counts in 
the UV and optical. Mazzei et al. (1992) were the 
first to try to model the evolution of stellar emis- 
sion and dust emission together in a consistent 



framework based on stellar population synthesis 
models and a physical calculation of dust absorp- 
tion. This model was then used by Franceschini et 
al. (1994) to calculate galaxy evolution and num- 
ber counts in bands from the optical through to 
the far-IR, based on the backwards evolution ap- 
proach. However, these models still made a num- 
ber of simplifying assumptions (e.g. slab geometry 
for disks), and set a number of present-day proper- 
ties of galaxies from observations (e.g. the optical 
depth of galactic disks, and the intensity of the 
radiation field heating the dust), rather than pre- 
dicting them. Simpler backwards evolution mod- 
els, where the luminosity evolution is parameter- 
ized as a simple function of redshift, have been 
considered by, e.g. Pearson & Rowan-Robinson 
(1996). 

Recently, dust absorption has been included 
in several different semi-analytical models (Kauff- 
mann et al. 1999; Somerville & Primack 1999; 
Cole et al. 2000). The first two calculate dust 
effects only for present-day galaxies, using a ID 
slab model, and taking the dust optical depth 
from observational measurements. On the other 
hand, Cole et al. (2000) predict the dust optical 
depth and how it evolves, based on chemical evo- 
lution and a prediction of disk sizes, and use the 
3D disk+bulge radiative transfer models of Fer- 
rara et al. (1999) to calculation the dust atten- 
uation. The only previous semi-analytical model 
to calculate dust emission as well as absorption is 
that of Guiderdoni et al. (1998). However, that 
model also has several limitations: the galaxy for- 
mation model does not include merging of either 
dark halos or visible galaxies, and the fraction of 
star formation occuring in bursts is simply an arbi- 
trary function; dust absorption is again modelled 
assuming a ID slab geometry; and the dust tem- 
perature distribution is not predicted. Instead, the 
dust emission spectrum is modelled as the sum of 
several components, whose temperatures and rel- 
ative strengths are chosen so as to reproduce the 
observed correlations of IR colours with IR lumi- 
nosity found by IRAS. 

The present paper represents a major advance 
over this earlier work in terms of scope, physi- 
cal self-consistency and predictiveness. We com- 
bine the semi-analytical galaxy formation model 
of Cole et al. (2000) with the stellar population 
+ dust model of Silva et al. (1998). The galaxy 
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formation model includes formation of dark halos 
through merging, cooling and collapse of gas in 
halos to form disks, star formation in disks regu- 
lated by energy input from supernovae, merging of 
disk galaxies to form elliptical galaxies and bulges, 
bursts of star formation triggered by these merg- 
ers, predictions of the radii of disks and spheroids, 
and chemical enrichment of the stars and gas. The 
stellar population + dust model includes a realis- 
tic 3D geometry, with a disk and bulge, two phase 
dust in clouds and in the diffuse ISM, star forma- 
tion in the clouds, radiative transfer of starlight 
through the dust distribution, a realistic dust 
grain model including PAHs and quantum heat- 
ing of small grains, and a direct prediction of the 
dust temperature distribution at each point in the 
galaxy based on a calculation of dust heating and 
cooling. The output is the luminosity and spec- 
trum of the stellar populations attenuated by dust, 
and of the dust emission from grains at a range 
of temperatures. From this, we can calculate the 
distribution of galaxy properties at any redshift, 
including the complete spectrum of each galaxy in 
the model from the far-UV to the sub-mm. 

In this paper we compare the predicted proper- 
ties for local galaxies with a wide range of obser- 
vational data. A future paper will be devoted to 
high-z galaxies (Lacey et al. 2000). In Sections 2 
and 3 we describe, respectively, the galaxy for- 
mation model and the spectrophotometric model. 
Section 4 describes how we generate model galaxy 
catalogues for both normal and starburst galaxies. 
The comparison with observations (SEDs, extinc- 
tion properties, colors, etc.) is presented in Sec- 
tion 5 for spiral galaxies, and in Section 6 for star- 
bursts. The model luminosity functions at differ- 
ent wavelengths are compared with observations 
in Section 7. Section 8 uses the models to predict 
the relationship between the star formation rate 
and the luminosities in various UV and IR bands, 
and to assess the accuracies of these as star for- 
mation indicators. Section 9 presents a summary 
and conclusions. 

2. Semi-analytical galaxy formation model 

We calculate the formation histories and global 
properties of galaxies using the semi-analytical 
galaxy formation model (GALFORM) of Cole et 
al. (2000) , a development of that described in Cole 



Table 1: Adopted values for GALFORM parame- 
ters. 



Cosmology 



o 


0.3 


A 


0.7 


h 


0.7 




0.02 


r 


0.19 


0"8 


0.93 



Star formation and feedback 



'.disk " 6.7 x 1(T 3 
a* -1.5 
ahot 2.0 
V bat (km/s) 150.0 

Stellar populations 
IMF Kcnnicutt (1983) 

T 1.4 
p 0.02 
R 0.29 

Mergers and bursts 

/ellip 0.3 
£*burst 0.5 
Tj 0.1 



et al. (1994) and Baugh et al. (1998). The princi- 
ple of the model is to calculate the formation and 
evolution of dark matter halos starting from an as- 
sumed cosmology and initial spectrum of density 
fluctuations, and then to calculate the evolution of 
the baryons (gas and stars) within these evolving 
halos using a set of simple, physically-motivated 
rules to model gas cooling, star formation, super- 
nova feedback and galaxy mergers. We describe 
here only the main features of the model, and re- 
fer the reader to Cole et al. (2000) for more details 
and for a discussion of the effects of varying pa- 
rameters with respect to standard values given in 
Table 1. 

(a) Cosmology: The cosmology is specified by 
the present-day density parameter fi , cosmolog- 
ical constant Ao, baryon fraction Qb (all in units 
of the critical density) and the Hubble constant 
Ho = lOOh kms~ 1 Mpc~ 1 . We assume a cold dark 
matter (CDM) model, with the initial spectrum 
of density fluctuations having shape parameter T 
and amplitude cr 8 (the r.m.s. density fluctuation 
in a sphere of radius 8/i _1 Mpc). 
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(b) Halo evolution: Dark matter halos form 
through a process of hierarchical clustering, build- 
ing up through merging from smaller objects. 
At any cosmic epoch, we calculate the number 
density of halos as a function of mass from the 
Press-Schechter (1974) formula. We then calcu- 
late halo merger histories, describing how a halo 
has formed, for a set of halos of different masses, 
using a Monte-Carlo algorithm based on the ex- 
tended Press-Schechter formalism. We generate 
many different realizations of the merger history 
for each halo mass. We then follow the process 
of galaxy formation separately for each of these 
realizations. 

(c) Cooling and collapse of gas in halos: 

Diffuse gas is assumed to be shock-heated to the 
virial temperature of the halo when it collapses, 
and to then cool radiatively out to a radius de- 
termined by the density profile of the gas and the 
halo lifetime. The gas which cools collapses to 
form a rotationally supported disk, for which the 
half-mass radius rdisk is calculated assuming that 
dark matter and associated gas are spun up by 
tidal torques, and that the angular momentum of 
the gas is conserved during the collapse. The gas 
supply by cooling is assumed to be continuous over 
the lifetime of the halo. 

(d) Star formation in disks: Stars form from 
the cold gas in the disk, at a rate 

ip = M co l d /T* disk , (1) 

where the star formation timescale is assumed to 
be 

T*disk = <Q 1 isk r disk (Vdisk/200kms _1 ) a * (2) 

where Vdisk is the circular velocity at the half- mass 
radius of the disk, and Tdisk = fdisk/ Vdisk is the dy- 
namical time. e*disk is the fraction of gas converted 
into stars in one dynamical time, for a galaxy with 
circular velocity Vdisk = 200km s -1 . The scaling of 
the star formation timescale with dynamical time 
is motivated by observations of star formation in 
nearby galaxies (Kennicutt 1998), but modified to 
reproduce the observed dependence of gas fraction 
on luminosity. 



(e) Supernova feedback in disks: The energy 
input from supernovae is assumed to reheat gas in 
the disk and eject it into the halo at a rate 

Mcjoct = /?diskV', (3) 

where for /?disk we assume 

/3disk = (V disk /V hot )- ahot (4) 

Gas which has been ejected is assumed to be un- 
available for cooling until the halo has doubled in 
mass through merging. The motivation for the 
this parameterization is that the rate of gas ejec- 
tion should be proportional to the rate of super- 
novae, and also depend on the escape velocity from 
the disk, which in turn is related to the circular 
velocity. Our standard case cthot = 2 is equiva- 
lent to the assumption that a constant fraction of 
the Type II supernova energy goes into ejecting 
gas from the disk, if the escape velocity is propor- 
tional to Vdisk- 

(f) Galaxy mergers and morphology: The 

galaxy morphology (i.e. whether it is a spiral or 
elliptical) is determined by merging. Following 
the merger of two halos, the largest pre-existing 
galaxy is assumed to become the central galaxy 
in the new halo, while the other galaxies become 
satellite galaxies. The central galaxy can continue 
to grow a disk by cooling of gas from the halo. The 
satellite galaxies merge with the central galaxy on 
a timescale equal to that for dynamical friction to 
make the orbits decay. The merger is classed as 
a major merger if the mass ratio of the satellite 
to central galaxy exceeds a value / e iii P , and as a 
minor merger otherwise. In a major merger, any 
pre-existing stellar disks are destroyed, producing 
a stellar spheroid (elliptical galaxy or bulge), and 
any remaining cold gas is consumed in a burst 
of star formation. The star formation timescale 
in the burst is related to the dynamical time of 
the bulge as described below. The spheroid can 
grow a new disk by cooling of halo gas. In a mi- 
nor merger, the stars from the satellite galaxy add 
to the bulge of the central galaxy, while the cold 
gas adds to the disk, but no burst is triggered. 
In either case, the half-mass radius rb u i g e of the 
spheroid produced in a merger is calculated using 
an energy conservation argument. Galaxies are 
classified into different morphological types based 
on their bulge-to-disk luminosity ratios. 
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(g) Star formation and feedback during 
bursts: As already mentioned, star formation 
bursts are assumed to be triggered by major merg- 
ers of galaxies. In Cole et al. (2000), these bursts 
were modelled in a very simple way, with the 
conversion of gas into stars being assumed to be 
instantaneous, since the galaxy properties exam- 
ined there were not sensitive to the detailed time 
dependence. In this paper, we model the bursts 
in more detail. We assume that star formation 
during bursts follows a law analogous to that for 
star formation in disks: 

1p = M C old/r*burst, (5) 

with star formation timescale 

Ttburst = e ^burst T bulgo (6) 

where Tbuige = ^buige/H>uige is the dynamical time 
of the spheroid formed in the merger, Vbuigc be- 
ing the circular velocity at rb u i go . As in Cole et 
al. (2000), feedback is modelled as in disks except 
with Vbuigc replacing Vdisk in eqn.(4), assuming the 
same values for Vhot and a no t, giving a feedback 
factor /3burst- Since we assume that no new gas is 
supplied by cooling during the burst, the star for- 
mation rate and cold gas mass decay during the 
burst as exp(— t/r e ), where 

T e = T*burst/ (1 - R + A>urst), (7) 

and R is the recycled fraction, discussed below. 
The burst is assumed to occur in a region of half- 
mass radius fburst, where 

n>urst = V r bu lgc (8) 

More details on the geometry assumed for star- 
bursts are given in § 3. For simplicity, the metallic- 
ity of the gas in the burst and of the stars formed 
during the burst are taken to be constant, and 
equal to the mean metallicity of the stars formed 
during the burst as calculated by the GALFORM 
model. The star formation in a burst is truncated 
at a time 5r e after the burst began, i.e. after 99% 
of the gas in the burst has either been converted 
into stars or blown out of the galaxy by super- 
nova feedback. At this time, the remaining gas 
and dust in the burst region are assumed to be 
dispersed. Star formation then starts again in a 
normal galactic disk surrounding the bulge, if one 
has formed by cooling of halo gas since the major 
merger that triggered the burst. 



(h) Chemical evolution: We assume that a 
fraction 1/T of the mass formed into stars goes 
into visible stars (0.1 < m < 125M ), while the 
remainder goes into brown dwarfs (m < O.IMq). 
For visible stars we adopt a universal IMF, similar 
to that in the solar neighbourhood. In Cole et al. 
(2000) and in this paper, we use the form proposed 
by Kennicutt (1983), which is consistent with the 
"best estimate" of Scalo (1998): 

dN/dlnm cx m~ ' 4 (m < 1M Q ) 

cx m- 1 - 5 (m > IMq) (9) 

We use the instantaneous recycling approxima- 
tion to calculate the evolution of the abundance 
of heavy elements of the cold gas (Z C oid) and stars 
(Z*) in each galaxy, together with that of the hot 
gas in the halo (Z hot ), including the effects of in- 
flows and outflows between the galaxy and halo. 
The chemical evolution depends on the recycled 
fraction R and the yield of heavy elements p. 

(i) Stellar population synthesis and dust ex- 
tinction: In Cole et al. (2000), we calculated 
the luminosity evolution of each galaxy at differ- 
ent wavelengths using the stellar population syn- 
thesis models of Bruzual & Chariot (2000). The 
effects of dust extinction were calculated in a sim- 
ple way using the dust models of Ferrara et al. 
(1999), which assume a smooth (unclumped) dis- 
tribution for both the dust (in a disk) and stars 
(in a disk and a bulge). In the present paper, we 
use instead the combined stellar population and 
dust model GRASIL (Silva et al. 1998) to calcu- 
late the galaxy luminosities and spectra including 
both extinction and emission by dust. The stel- 
lar population part of GRASIL is similar to the 
Bruzual & Chariot model, as both are based on 
similar stellar evolution tracks and stellar spectra. 
The dust part of GRASIL is however considerably 
more sophisticated than the Ferrara et al. mod- 
els, in that GRASIL allows for clumping of both 
dust and stars, and calculates the grain heating 
and emission as well as the extinction. 

The parameters we have chosen for the GAL- 
FORM model are the same as those of the stan- 
dard AC DM model of Cole et al. (2000), apart 
from e*burst and r\ describing the timescale and ra- 
dius of bursts, which were not considered in Cole 
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et al. . These parameters values are given in Ta- 
ble 1, and were obtained by comparing the model 
to observations of nearby galaxies, without any 
consideration of the far-IR or UV properties. We 
refer the reader to Cole et al. for a complete discus- 
sion of the effects of varying the 'old' GALFORM 
paramaters, and for a systematic presentation of 
the influence of these parameters on the optical- 
NIR properties of galaxies (LFs, Tully-Fisher re- 
lation, disk sizes, morphology, gas content, metal- 
licity, M/L ratios and colours). We only recall 
here the main observational constraints used to 
fix each of these parameters: e*disk - gas fraction 
of L± galaxies; a* - variation of gas fraction with 
luminosity; aihot - faint end of LF and Tully-Fisher 
relation ; Vh t - faint end of LF and sizes of low- 
ly spirals; IMF - observations of solar neighbour- 
hood; T - in LF; p - metallicity of L* ellipticals; 
/eiiip - morphological mix of galaxies. 

Values for e*burst and 77 are instead obtained 
later in this paper by detailed comparison of the 
results of the combined GALFORM+GRASIL 
models with observed properties of bursting galax- 
ies. 

3. The Stellar Population and Dust Model 

Far-UV to mm SEDs of model galaxies are 
calculated using the GRASIL code (Silva et al. 
1998), which follows both the evolution of stel- 
lar populations and absorption and emission by 
dust. GRASIL calculates the following: (i) emis- 
sion from stellar populations; (ii) radiative trans- 
fer of starlight through the dust distribution; (hi) 
heating and thermal equilibrium of dust grains (or 
thermal fluctuations for small ones); and (iv) emis- 
sion by dust grains. 

3.1. Stellar Population Model 

The single stellar population (SSP) libraries in- 
cluded in GRASIL are based on the Padova stellar 
models and cover a large range in age and metal- 
licity. They include the effects of dusty envelopes 
around AGB stars (Bressan et al. 1998). The age 
and metallicity distribution of a composite stellar 
population is specified by the birthrate function 
^(t, Z), where ^(t, Z) dt dZ gives the mass of stars 
that were formed in the time interval (t, t+dt) with 
metallicities in the range (Z, Z + dZ) . The SED 
for the composite stellar population at time t is 



then obtained using 



Lx(t) 



dt' / dZ l x (t 
Jo 



t',Z)H>(t',Z) (10) 



where 1\(t, Z) is the SED of a SSP of age t and 
metallicity Z for the assumed IMF. 

For our semi-analytical galaxy formation model, 
^(t, Z) is calculated for each galaxy by summing 
over all the progenitor galaxies which have merged 
to produce that galaxy, separately for the disk and 
bulge components. The progenitor galaxies each 
had their own star formation and chemical history, 
so that the composite ^(t, Z) obtained in general 
has a broad distribution of metallicity at each each 
age, i.e. there is no unique age-metallicity relation 
Z(t). 




diffuse ISM 
+ free stars 



BULGE 



Fig. 1. — Sketch of geometry of stars and dust in 
the GRASIL model. 



3.2. Dust Model 

GRASIL computes the radiative transfer of 
starlight, the heating of dust grains, and the emis- 
sion from these grains with a self-consistent calcu- 
lation of the distribution of grain temperatures, for 
an assumed geometrical distribution of the stars 
and dust, and a specific grain model. The dust 
is divided into two components, dense molecular 
clouds and diffuse cirrus in a disk. Stars form in- 
side clouds and progressively leak out. 
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Table 2: Adopted values for the adjustable 
GRASIL parameters. Note that the results de- 
pend on the ratio M c /r 2 rather than on the two 
quantities individually (sec § 3.2). 



h z /h R {disk) 


0.1 


h z /hn{burst) 


0.5 


h z {dust)/ h z (stars) 


1 


fmc 


0.25 


Mjrl 


10 6 M Q /(16pc) 2 


tesc(disk) 


2Myr 


t esc (burst) 


lOMyr 



The details are given in Silva et al. (1998), 
but for convenience we summarize the main fea- 
tures here, focusing on the modifications intro- 
duced for the purposes of this application. Those 
GRASIL parameters which are not provided by 
GALFORM, and are in this sense additional 
adjustable parameters of the combined GAL- 
FORM+GRASIL semi-analytic modelling, are 
listed in Table 2, together with the adopted values 
for our standard case. See Fig. 1 for a sketch of 
the geometry of our model. 

(a) Geometry of stars: The stars are in two 
components (Silva et al. (1998) considered only 
pure disk and pure bulge systems): (i) a spher- 
ical bulge with an analytic King model profile, 
p oc (r 2 + r 2 ) -3 / 2 for r < r t , with concentra- 
tion parameter log(r t /r c ) = 2.2; (ii) a disk with 
a radially and vertically exponential profile, scale- 
length hn and scaleheight h z . As described in § 2, 
the disk and bulge masses, M<ji s k and Mbuige; and 
half-mass radii, rdi s k and rtmige, for any galaxy 
are predicted by the galaxy formation model. The 
bulge core radius is related to the bulge half-mass 
radius by r c = rb u i g c/14.6, while the disk scale- 
length Jir is related to the disk half-mass radius by 
hit = fdisk/1-68. The star formation histories are 
also calculated separately for the disk and bulge by 
GALFORM. However, the disk axial ratio h z /hn 
is a free parameter of the GRASIL model. 

As partially anticipated in § 2, in galaxies un- 
dergoing bursts, the burst star formation, as well 
as the gas and dust, are assumed to be in an ex- 
ponential disk, but with half-mass radius rburst = 
V r bui g o rather than rdisk- The axial ratio h z /hu 
of the burst region is allowed to be different from 



that for disks in non-bursting galaxies. The stars 
which were formerly in the disks of the galaxies be- 
fore the galaxy merger which triggered the burst 
are assumed to become part of the bulge following 
the merger. 

(b) Geometry of gas and dust: The gas and 

dust are in an exponential disk, with the same ra- 
dial scalelength as the disk stars (either rdi s k for 
normal galaxies or rburst = V n>uige for starbursts), 
but in general with a different scaleheight, so that 
h z (dust) /h z (stars) is a free parameter. The gas 
and dust are in two components within the disk, 
molecular clouds and the diffuse ISM. The latter 
corresponds to the cirrus dust. The total gas mass 
Mcoid and its metallicity Z C oid are calculated by 
the galaxy formation model, but the fraction of 
the gas in clouds, f mc , and the cloud mass M c and 
radius r c are free parameters of GRASIL, though 
the results actually depend only on their combina- 
tion M c /r 2 , which determines, together with the 
dust/gas ratio (see point (d) below), the optical 
depth of the clouds (Silva et al. 1998). 

(c) Young stars and molecular clouds: Stars 
are assumed to form inside the molecular clouds, 
and then to escape on a timescale t csc . Specifically, 
the fraction of stars still inside clouds at time t 
after they formed is assumed to be given by 

F(t) = 1 (t<t csc ) 

= (t>2t csc ) (11) 

We allow t csc to take different values in normal 
disks and in bursts, in keeping with the results 
of Silva et al. (1998). Indeed, given the small 
size scale and the intensity of the star formation 
activity in bursts, it is conceivable that the star- 
forming environment is quite different from that 
in normal spiral galaxies (see also § 3.3). 

(d) Dust grain model: The dust is assumed to 
consist of a mixture of graphite and silicate grains 
and polycyclic aromatic hydrocarbon molecules 
(PAHs), each with a distribution of grain sizes. 
Absorption and emission properties are calculated 
for each grain composition and size. The grain 
mix and size distribution were chosen by Silva 
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et al. (1998) to match the extinction and emis- 
sivity properties of the local ISM, and arc not 
varied here. The dust/gas ratio 5 in the clouds 
and diffuse ISM is assumed to be proportional to 
the gas metallicity, with a value S — 1/110 for 
Z = Zq = 0.02. Thus, the total dust mass in a 
galaxy scales as M dust oc Z co \ A M cold . 

(e) Radiative transfer, dust heating and 
re-emission: The luminosities of the different 
stellar components (bulge stars, disk stars, and 
young stars still in clouds) are calculated using 
the population synthesis model described above. 
The GRASIL code then calculates the radiative 
transfer of the starlight through the dust distri- 
bution. Whilst in molecular clouds a full radia- 
tive transfer calculation is performed, the effects 
of scattering by diffuse dust are included only ap- 
proximately, by assuming that the effective optical 
depth for absorption is related to the true absorp- 
tion and scattering optical depths T a b s and T sca t 
by T abSie ff = ^ T abs (T abs + T scat ) . Thus the dust- 
attenuated stellar radiation field can be calculated 
at any point inside or outside the galaxy. Then 
GRASIL calculates for each point in the galaxy 
the absorption of radiation, thermal balance and 
re-emission for each grain composition and size. 
Thus, the distribution of grain temperatures is cal- 
culated self-consistently for the assumed geometry 
of the stars and dust, including the effects of tem- 
perature fluctuations for small grains. The final 
galaxy SED L\ is obtained by adding the con- 
tributions from the starlight (attenuated by dust) 
and from the dust re-emission, and depends on the 
inclination angle at which the galaxy is viewed. 
Emission from dust in the envelopes of AGB stars 
is included in the SSPs. 

Our computations allow us to calculate the 
amount of energy emitted in the PAH bands, but 
theoretical predictions of the detailed shapes of the 
emission features are rather uncertain. Therefore 
we use the Lorentzian analytical fits to the ob- 
served PAH profiles for the Ophiuchus molecular 
cloud from Boulanger et al. (1998). 



3.3. Choice of GRASIL adjustable param- 
eters and new GALFORM parame- 
ters 

The values of GRASIL parameters (Table 2) 
not provided by GALFORM have been based on 
a variety of observational data for galaxies in the 
local universe. For some of them, the choices were 
made by trying to match model predictions to the 
observational data, as is discussed in more detail 
in the relevant sections of this paper. We now 
summarize the reasons for these choices and for 
those of the two GALFORM parameters (e*burst 
and 77) not considered in Cole et al. (2000) 

(a) e^burst^ this is chosen mainly so as to repro- 
duce the bright end of the IR luminosity function, 
which is dominated by bursts triggered by galaxy 
mergers (§7.4). A secondary (weak) constraint is 
to reproduce the relation between Lm/Luv an d 
total luminosity or UV slope j3 observed for star- 
burst nuclei (§6.2). The value controls both the 
luminosity and lifetime (and thus number density) 
of starbursts. 

(b) r] — rburst/^buigo: the choice of this is mainly 
based on the observational fact that starburst ac- 
tivity is usually confined to a nuclear region with 
a size much smaller than the galaxy as a whole, 
by about one order of magnitude (e.g. Sanders 
& Mirabel 1996, and references therein). For 
instance, in Arp 220 most of the molecular gas 
is found in the central ~ 300pc (Scoville et al. 
1997), and the mid-IR light is dominated by more 
or less the same region (Keto et al. 1992), while 
the half-light radius for the old stellar population 
is - 3kpc (Wright et al. 1990). The value of 
T] controls the amount of extinction of starlight 
from bursts by the diffuse ISM, which however is 
usually overwhelmed (in bursts) by extinction in 
molecular clouds (see §6.4). Therefore our results 
are not very sensitive to the precise choice of this 
parameter, nor to the value of h z /h R in starbursts 
(discussed below). 

(c) h z /IiR: for normal disks, we choose a value 
of 0.1, consistent with observations of the stel- 
lar light distributions in edge-on spiral galaxies 
(e.g. Xilouris et al. 1999). It is also the typ- 
ical value used by Silva et al. (1998) to fit the 
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SEDs of spiral galaxies. This value is also im- 
portant, and the adopted value turns out to be 
suitable, to match the observed difference in ex- 
tinction between spiral galaxies seen edge-on and 
face-on (§5.2). Apart from this test, most pre- 
dicted properties are not very sensitive to h z /hn- 
The choice of h z /hn = 0.5 for starbursts is based 
on general observational indications that they are 
only moderately flattened. 

(d) h z (dust)/h z (stars): this parameter has a 
significant effect on how much starlight is ab- 
sorbed in the diffuse medium. From observations 
of our own galaxy it is known that the scale- 
height of stars increases with the age of the stellar 
population, so that there is no unique value for 
h z (dust)/h z (stars). The scaleheight of the gas is 
comparable to that of the youngest stars. Since 
we are particularly interested in having a real- 
istic estimate of the extinction in the UV, both 
because it is strongest there and because this is 
an important source for dust heating, we choose 
h z (dust) I h z (stars) = 1 to match what is seen for 
the young stars. 

(e) f mc : this can be estimated observationally 
from the ratio of molecular to atomic hydrogen in 
galaxies, since in normal spiral galaxies most of 
the hydrogen in molecular clouds is in H2, while 
most of the intercloud medium is atomic HI. Our 
adopted f mc implies a ratio H 2 /HI similar to the 
typical one for L* spirals found by Sage (1993). 
Larger values reduce the extinction in the diffuse 
ISM and produce a somewhat colder molecular 
clouds emission, but our results are in general not 
significantly affected as long as we keep f mc in the 
range 0.2-0.8. 

(f ) M c ,r c : as already remarked (§ 3.2) the pre- 
dicted SEDs depend on the ratio M c /r^, rather 
than M c and r c separately. Thus M c has been 
chosen to match typical giant molecular clouds in 
our own and nearby galaxies, while r c is chosen 
based on the results of Silva et al. (1998), who 
tuned M c /rj? to fit the SEDs of starburst galaxies 
in particular. The resulting value for r c is consis- 
tent with direct measurements of cloud radii. 



(g) ^csc: this is a very important parameter in 
the model, since it is this that mainly controls how 
much of the radiation from young stellar popula- 
tions is absorbed by dust. Silva et al. (1998) found 
from detailed fits to 3 nearby spirals values of 2.5, 
3 and 8 Myr. For normal spirals, we favor a value 
of 2 Myr, close to the lower limit of this range, 
rather than the average 5 Myr. Although the lat- 
ter provides an equally good overall fit to the LFs 
(somewhat better for IRAS colors and LFs, § 7.4, 
but somewhat worse for the UV LF), the former 
is more consistent with the massive star census 
in our own and nearby galaxies, which suggests 
that the time for which the stars are obscured by 
dust is about the 20 % of the total lifetime for 
the brightest stars, above say 30 M©, whose life- 
time is around 6 Myr. For starbursts, the value 
we choose is based mainly on the comparison with 
properties of UV-bright starbursts in §6.2. This 
leads us to a value closer to that of normal spirals 
than the values t esc = 20-60 Myr found by Silva 
et al. (1998) from fitting 3 nearby starbursts, and 
suggests that the starburst galaxies used by Silva 
et al. may be not representative of the whole popu- 
lation. The difference could also be due in part to 
the more complex geometry adopted in this paper 
for starburst galaxies. 

4. Generation of model galaxy catalogues 

The GALFORM model is run for a set of dark 
matter halos covering a large range in mass, and 
generates a catalogue of model galaxies, includ- 
ing information about the following properties for 
each galaxy at the chosen epoch: stellar masses 
M disk , Mbuigo, and half mass radii r disk and r bu i g o 
of the disk and bulge, mass M co id and metallicity 
■Zcoid of gas in the disk, and the star formation 
histories ^(t, Z) of the disk and bulge separately, 
including both star formation in disks and during 
bursts, and specifying the metallicity distribution 
of the stars of each age. In addition, each galaxy 
has a weight or number density n, such that that 
galaxy should appear N — nV times in an average 
volume of the universe V. 
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The GALFORM code outputs all the galaxies 
for each different halo that is calculated, down to 
a minimum mass controlled by the mass resolu- 
tion of the merger tree. In practice, this means 
that the model catalogue contains many more low 
mass galaxies than high mass galaxies. Running 
the GRASIL code on every galaxy in the original 
catalogue is neither feasible (because of computer 
time) nor necessary. We therefore select a subset 
of galaxies from the catalogue chosen to sample 
galaxies more evenly in mass, and redistribute the 
weights to give the same total number density in 
each mass range. The GRASIL code is then run on 
each galaxy in this reduced catalogue to give the 
SED L\ including both stellar emission and dust 
absorption and emission, and statistical proper- 
ties (e.g. luminosity functions) are then calculated 
making use of these weights. In fact, we calculate 
2 samples of galaxies, a "normal" sample and a 
"burst" sample, as follows: 

(a) Normal galaxies: By "normal" galaxies, we 
here simply mean galaxies not selected to have had 
a recent burst. From the parent GALFORM cat- 
alogue, we select a sample with equal numbers of 
galaxies in equal bins in log M„ , M* being the to- 
tal stellar mass of the galaxy. Within each mass 
bin, galaxies are randomly selected (allowing for 
multiple selection of the same galaxy) with prob- 
ability proportional to their weight n. The se- 
lected galaxies are then assigned new weights rij, 
such that each galaxy within the same bin has 
the same weight (multiply selected galaxies being 
counted as separate objects), and that the sum of 
the weights (i.e. number densities) within a bin is 
the same as in the parent catalogue. We have used 
bins with AlogM, = 0.3 and about 40 galaxies 
per bin. 

(b) Burst galaxies: By "burst" galaxies we mean 
galaxies which have had a burst in the recent 
past, at whatever redshift we are looking. Bursts 
have short durations compared to the age of the 
universe, so the fraction of galaxies undergoing a 
burst at any one time is very small, but they can 
be very luminous, and so may dominate the galaxy 
luminosity function at the highest luminosities. In 
practice, our "normal galaxy" catalogue contains 
too few galaxies in total to provide a represen- 
tative sample of galaxies seen during their burst 
phase. Rather than use a greatly enlarged "nor- 
mal galaxy" sample, it is more efficient to calcu- 



late a separate sample of "burst" galaxies, as fol- 
lows: for a redshift z, we choose a subsample of 
galaxies which have had bursts during the time in- 
terval t(z) > t > t(z) — T, where t(z) is the age 
of the universe at redshift z, with equal numbers 
of galaxies in equal bins in logMburst, -Mburst be- 
ing the mass of stars formed in the most recent 
burst. The galaxies are assigned new weights rij 
analogously to the case of normal galaxies, but 
now conserving the total number density in bins 
of Mburst for the galaxies which have had bursts 
more recently than T. For each burst galaxy, we 
then run GRASIL to calculate the total galaxy lu- 
minosity at a set of times after the start of the 
burst, chosen to sample all phases of the burst 
evolution, including the highest luminosity phase 
of short duration. If T « t(z), then the rate of 
bursts per unit volume during the time interval T 
can be taken as constant. Then, for the ith galaxy 
in the jth phase in the burst evolution that lasts 
a time Atj, the number density of galaxies that 
should be found in this phase is 




(12) 



These weights can then be used to calculate sta- 
tistical properties such as luminosity functions. 
When combining the "normal" and "burst" galaxy 
samples, the normal galaxies with bursts more 
recent than T are explicitly excluded, to avoid 
statistical double-counting. In practice, we chose 
T = t{z)/2Q at all z, with bins AlogM bu rst = 0.3, 
around 10 galaxies per bin, and around 10 out- 
put times per galaxy, for < t — hurst < 100r e . 
For many calculations of statistical distributions, 
we then interpolate between these output times to 
have more burst phases. 

5. Properties of spiral galaxies 

In this section, we test the model predictions 
for disk galaxies against observed emission and ab- 
sorption properties of nearby spirals. 

5.1. SEDs of face-on spirals 

We compared the predicted near-UV to far-IR 
SEDs of our model galaxies with the broad-band 
SEDs of a complete sample of nearby spiral galax- 
ies (de Jong & van der Kruit 1994), consisting 
of a diameter-limited sample of 86 nearly face-on, 
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Fig. 2. — Upper panels: SEDs (normalized to the K-band) of a representative sample of spiral galaxies from 
the models. Middle panels: Comparison of the models with the observed SEDs of a sample of face-on spiral 
galaxies obtained by de Jong & van der Kruit (1994). The median and the 10% and 90% percentiles of the 
model SEDs are plotted as lines. Lower panel: comparisons between observed SEDs of individual galaxies 
from the sample and selected model SEDs. The left-hand panels are for galaxies in the luminosity range 
-23 > M K - 51og/i > -24, and the right-hand panels for -22 > M K - 51og/i > -23. 
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disk-dominated galaxies, de Jong & van der Kruit 
measured fluxes of these galaxies in the BVRIHK 
bands. We have supplemented these with U- 
band magnitudes from the literature and IRAS 
12, 25, 60, 100/rni fluxes from Saunders (1997). 

We considered only those model galaxies with 
bulge to total light ratio B/T < 0.5 in the B-band, 
corresponding to the range of types in the de Jong 
sample (e.g. Simien & de Vaucouleurs 1986). 
From Figure 2 it is apparent that the models re- 
produce the observed spectral trends reasonably 
well. This is particularly impressive since the ra- 
tio between the infrared and the optical-UV spans 
more than one order of magnitude, both in the ob- 
served and in the theoretical SEDs. The predicted 
infrared emission peaks at wavelengths somewhat 
larger than those sampled by IRAS, in agreement 
with recent ISO observations (e.g. Alton et al. 
1998). The emission in the mid- infrared is domi- 
nated by PAH molecular bands. 

Figure 3 shows the effects on typical SEDs of 
factor 2 variations in the molecular cloud fraction 
f mc , their mass M c (i.e. their optical depth, hav- 
ing fixed the radius) and the escape timescale t csc . 
The effects are mostly confined to the mid-IR be- 
tween 8 and 40 /jm and in the UV below 0.4 /xm. 
In these spectral regions, the predicted flux may 
change by up to a factor ~ 2, while the effects 
are almost negligible elsewhere. In the mid-IR 
the most important parameter is the cloud optical 
depth, while in the UV the effects of t osc dominate. 
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Fig. 3. — Effects on typical SEDs of factor 2 vari- 
ations of f mc , M c (effectively, the cloud optical 
depth) and t csc . In each case only one parameter 
is varied with respect to the standard model. 




M K — 51ogh 

Fig. 4. — Predicted dependence of dust opti- 
cal depth on luminosity for spiral galaxies in the 
model, tvo is the central face-on extinction optical 
depth in the V-band for the diffuse dust compo- 
nent in the disk. The line shows the median, and 
the error bars show the 10% and 90% percentiles. 
Model galaxies are selected to have B/T < 0.4 in 
the B-band. 
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5.2. The global extinction in spiral galax- 
ies 

The models predict that the extinction in 
galaxy disks should increase strongly with galaxy 
luminosity as shown in Figure 4. Clearly in com- 
paring predictions of dust extinction with observa- 
tions, one must be careful to specify the luminosity 
of the objects concerned. 

There have been many attempts to measure or 
observationally constrain the total dust extinction 
in galaxy disks, using a variety of techniques: the 
inclination dependence of magnitudes or colours 
(e.g. de Vaucouleurs et al. 1991; Giovanelli et al. 

1995) , surface brightness distributions in edge-on 
galaxies (e.g. Kylafis & Bahcall 1987), colour 
gradients in face-on disks (e.g. de Jong 1996c), 
and the ratio of far-IR to UV luminosities (e.g. Xu 
&Buat 1995;Buat&Xu 1996; Wang & Heckman 

1996) . In general, different techniques have given 
somewhat different answers. 

Xilouris et al. (1999) estimate dust extinctions 
by fitting detailed models of the star and dust dis- 
tributions to the observed surface brightness dis- 
tributions of edge-on spiral galaxies. Their dust 
models include scattering. From their six Sb-Sc 
spirals with luminosities in the range —17.5 > 
Mb — hlogh > —19.0, we obtain a median cen- 
tral face-on extinction optical depth ryo = 0.6. In 
comparison, for edge-on galaxies in our model in 
the same luminosity range (after extinction), and 
with 0.1 < B/T < 0.3 in the B-band, we find 
a median value Ty = 2.2, which is significantly 
larger. There could be several reasons for this dif- 
ference between the models and the observations: 
there may be a problem with the Xilouris et al. 
method for deriving ryo from the observations, or 
the Xilouris et al. sample may not be representa- 
tive, or the problem might be with our assumption 
that the dust and stars have the same exponen- 
tial scalelength. The extinction-inclination obser- 
vational test discussed next implies extinctions for 
edge-on galaxies in this luminosity range which are 
at least as large as those predicted by our model. 

We considered the dependence of the net ex- 
tinction on the inclination angle at which a galaxy 
is viewed. This has been studied in many pa- 
pers using different methods, most recently by 
Tully et al. (1998), who also summarize the re- 
sults from the earlier studies. Tully et al. mea- 



sure the dependence of B — K, R — K and I — K 
colours on galaxy inclination at a given K-band 
luminosity the K-band being chosen to minimize 
extinction effects. They have a complete sam- 
ple of spirals covering a large range in luminos- 
ity -18.5 > M K - 5 log h > -24.5. They find 
a strong luminosity dependence, with a difference 
in B-band extinction between edge-on and face-on 
galaxies of about 2 mag for the brightest galaxies, 
and negligible for the faintest ones. 

Tully et al. (1998) follow the usual practice and 
parameterize the extinction relative to that for the 
galaxy seen face-on as 

4"° ee m x (i) - m A (0) = 7A log(a/6) (13) 

where 7a is a function of the passband. The axial 
ratio a/6 is assumed to be related to the inclination 
angle i by 

where i = for a face-on system, and q is the axial 
ratio of a galaxy seen edge-on. 

The models are compared with observations in 
Fig. 5. We use equation (14) to convert from the 
model inclination angle to the axial ratio, assum- 
ing q — 0.1, which is the ratio h z /fiR adopted in 
our galaxy models. We considered model galaxies 
corresponding to the morphological types Sa-Scd, 
and four ranges in K-band luminosity, correspond- 
ing to the ranges chosen by Tully et al. (1998), 
indicated by the different symbols in the figure. 

The three lines in the figure correspond to dif- 
ferent values of the slope 73. The model galax- 
ies approximately follow the linear dependence on 
log(a/6) (equation (13)), but with slopes jb that 
are somewhat shallower, at any given luminosity, 
than those observationally inferred by Tully et al. 
(1998). For instance, for the luminosity range 
-23.0 < M K - 51ogh < -22.0, our models fol- 
low an average slope (7s) ~ 0.9, while Tully et al. 
find jb = 1.1 ± 0.5, after allowing for the K-band 
extinction. The slope predicted by our models de- 
pends on the value chosen for the parameter h z /hn 
(see also § 3.3). We have checked that increasing 
h z /hft from our adopted value of 0.1 decreases the 
slope, while reducing it does not increase the slope 
significantly. 
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Fig. 5. — Dependence of extinction on inclination 
angle predicted by our galaxy models. The quan- 
tity A'g = M B (i) - M B (0) is plotted against the 
apparent axial ratio a/b, for the Sa-Scd models 
(with 0.04 < B/T < 0.4 in the B-band). The 
galaxies are plotted with different symbols in four 
K-band luminosity ranges, as detailed in the fig- 
ure. The same galaxies arc plotted seen at differ- 
ent inclinations. For clarity, small horizontal off- 
sets have been applied to the galaxies in different 
luminosity ranges. Three lines are plotted corre- 
sponding to equation (13) with the slopes 75 = 
0.5, 1 and 1.5. 



greater than its time-averaged value, and the star 
formation timescale correspondingly much shorter 
than the age of the universe. This definition in- 
cludes objects with a wide range of properties, 
from bursting dwarf irregular galaxies (e.g. Thuan 
& Martin 1981) to the ultra-luminous IR galaxies 
(ULIGs) found by IRAS (e.g. Sanders & Mirabel 
1996). In practice, a large variety of observational 
criteria have been applied to select samples of star- 
burst galaxies, ranging from optical morphologies 
and spectra (e.g. Balzano 1983) to IR colours and 
luminosities (e.g. Armus et al. 1990; Lehnert & 
Hcckman 1995). In our galaxy formation model, 
bursts are assumed to occur following major merg- 
ers of galaxies, producing elliptical galaxies from 
disk galaxies. For the ultra-luminous IR galax- 
ies, the link between the starburst activity and 
galaxy mergers is clearly established (e.g. Sanders 
& Mirabel 1996), while for low-luminosity star- 
bursts, additional triggering mechanisms probably 
operate, which arc not included in our model. In 
this section we will compare the properties of star- 
bursts predicted by our model with those of vari- 
ous observational samples. 



The agreement for the inclination test is any- 
way acceptable. Part of the discrepancy could be 
due to our simplified treatment of scattering by 
the diffuse dust (§ 3.2). In our models, the abso- 
lute extinction is often dominated by the molecu- 
lar clouds, but the difference between the face-on 
and the edge-on extinction is entirely due to the 
cirrus. Comparisons of our model with that of Fer- 
rara et al. (1999), where the treatment of scatter- 
ing is more accurate, show that for the brightest 
objects this effect can account for about 0.1-0.2 
mag of the differential extinction in the B band. 

6. Properties of starburst galaxies 

Starburst galaxies are broadly defined as galax- 
ies in which the current star formation rate is much 



15 



6.1. Properties of starbursts in the model 

Figure 6 shows how various properties of the 
bursts in our model vary with the stellar mass of 
the galaxy, M sta r, after completion of the burst. 
The total mass of new stars formed in the burst, 
■^burst) is seen to increase with the galaxy mass, 
with the fraction of stars formed in the burst be- 
ing typically between ~ 1% and ~ 50%. An ex- 
ception to this trend is the group of points in 
the lower right corner of Fig. 6a corresponding to 
small bursts occuring in large galaxies. These 
small bursts are produced by mergers between 
gas-poor elliptical galaxies. The main trend in 
panel (a) is produced by mergers between disk 
galaxies containing significant fractions of gas, 
and these dominate the statistics at all burst 
masses. The star formation rate during the burst 
is (Mburst/T e ) exp(— t/r e ), with t measured from 
the start of the burst. The peak star formation 
rate is thus Mb urs t/r e , and occurs at the begin- 
ning of the burst. This peak SFR is seen also 
to increase with the host galaxy mass. The half- 
mass radius rburst and exponential decay time T e 
of the burst are assumed to scale with the half- 
mass radius and dynamical timescale of the host 
galaxy, and also increase with galaxy mass. Large 
bursts, with Mburst ~ W w h~ 1 M ( ? ) , are predicted 
to occur in galaxies with M star ~ 10 11 M Q , and 
to have radii rburst ~ 0.5/i kpc, star formation 
timescales r e ~ 5 x 10 7 yr, and peak star forma- 
tion rates ~ 200/i~ 1 M Q yr~ 1 . These are similar 
properties to those inferred observationally for the 
ULIGS (e.g. Sanders & Mirabel 1996). 

6.2. Properties of UV-bright starbursts 

A large amount of work has been done on sam- 
ples of UV-bright starbursts selected from the cat- 
alogue of UV spectra of star-forming galaxies of 
Kinney ct al. (1993). Various correlations have 
been found, for instance between the bolomctric 
luminosity, the UV/IR ratio, the slope of the UV 
continuum and the mctallicity (e.g. Meurer et al. 
1995; Heckman et al. 1998). In this section, we 
compare the properties of our model starbursts 
with some of this observational data. 

The observational sample that we use for our 
comparison is that of Heckman et al. (1998), 
who selected 45 starburst and star-forming galax- 
ies from the original atlas of Kinney et al. . The 



criteria for a starburst galaxy to appear in the 
Kinney et al. catalogue are (a) that it has been 
previously classified as a starburst based on op- 
tical data, usually meaning that it has a com- 
pact optical morphology and strong optical emis- 
sion lines (but no AGN activity) (e.g. Balzano 
1983); and (b) that it has been observed by IUE 
and has a high enough surface brightness within 
the IUE aperture to produce a reasonable quality 
UV spectrum. The catalogue is not in any sense 
statistically complete. The starburst activity in 
these galaxies is generally confined to the central 
regions. (The galaxies have mostly been selected 
so that the starburst activity fits within the IUE 
aperture, 20" x 10", while the optical diameters 
of the underlying galaxies are typically a few ar- 
cminutes.) 

For the galaxies in their sample, Heckman et 
al. (1998) measured a UV luminosity Lyy = 
\L\(1900A) and mean continuum slope (3 between 
1250 and 1850 A (defined by L\ cx X 13 ) from IUE 
spectra, and a far-IR luminosity L fir from IRAS 
measurements. Heckman et al. use the definition 
of Lfir from Helou et al. (1988), which can be 
expressed as in terms of the luminosities in the 60 
and 100 /im IRAS bands as 

L FIR = 0.65^(60) + 0.421^(100) (15) 

Lfir provides an estimate of the 40 — 120/im 
luminosity. The quantity Lfir + Lyy is simi- 
lar to the bolometric luminosity in the case of 
starbursts, where most of the radiation is emit- 
ted in cither the UV or the FIR. The galaxies 
in the Heckman et al. sample cover the range 
Lfir + L uv ~ lO 8 -lO n L . 

The evolutionary tracks of a selection of 
model starbursts, with burst masses covering the 
range M burst ~ 10 7 - 1O 1O M , in Lfir + L uv , 
Lfir/Ljjv and (3 are shown in Figure 7, together 
with observational data for the Heckman et al. 
sample. We have calculated these quantities from 
the model SEDs to match the way they are cal- 
culated from the observational data. The bolo- 
metric luminosities of the model bursts, as mea- 
sured by Lfir + Luv, peak soon after the start 
of the burst, following which they evolve towards 
smaller values. At the same time, the amount of 
dust reprocessing of the radiation, as measured by 
Lfir/Luv, also decreases. This results from two 
effects: the escape of young stars from the dense 
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Fig. 6. — Properties of starbursts in the model at z = 0. (a) The mass of new stars formed in the burst, 
-^burst vs the stellar mass of the galaxy after the burst, M star . (b) The peak star formation rate in the burst 
-Mburst/T e vs M star . (c) The half-mass radius of the burst rburst vs M star . (d) The exponential decay time of 
the burst r e vs M star . 



molecular clouds, and the decrease in the optical 
depth of the diffuse dust component as the gas in 
the burst is consumed. The UV slope f3 initially 
evolves towards more negative values, i.e. bluer, 
as the net dust opacity falls. However, as the rate 
of formation of new stars drops and the dominant 
stellar population becomes older, the intrinsic un- 
absorbed stellar spectrum becomes redder, so the 



evolution in j3 reverses, the models becoming red- 
der with time even though the dust attenuation 
is falling. This happens after 20 — 30Myr, con- 
trolled mainly by the stellar evolution timescale. 
As long as the evolution in /3 is dominated by the 
declining dust opacity, the models stay close to 
the locus of observed points in the Lfir + Ljjy vs 
Lfir/Luv an d Lfir/Luv vs (3 panels, but when 
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Fig. 7. — Upper panel: Lpm + Ljjv (roughly 
corresponding to the bolometric luminosity) in 
/j-2 cr g S -i versus the ratio of the IR to the UV 
luminosity Lfir/ Luv ■ Lower panel: the rela- 
tion between Lfir/Luv and the slope (3 of the 
UV continuum. In each panel, the filled sym- 
bols show observed local starbursts from Heck- 
man et al. (1998), while the lines show the evo- 
lutionary tracks of five different model starbursts. 
The bursts start towards the upper right part of 
each panel, and evolve towards the lower left (top 
panel) or lower right (bottom panel). The crosses 
mark the age since the start of the burst in steps 
of 5 Myr. 

the intrinsic stellar spectrum starts to redden with 
age, the models move away from the observed lo- 
cus. This is not in itself in contradiction with 
observations, since there are selection effects in 
the observational sample, as we discuss below. 

The burst evolution involves the interplay be- 
tween two timescales, the lifetime of massive stars, 
~ 10 7 yr, and the exponential decay time r e of the 
star formation rate and gas mass in the burst. The 
latter varies with burst mass, being larger than the 
stellar evolution timescale for large bursts, and 
comparable for small bursts, as shown in Fig 6. 
The model starbursts begin their evolution with a 
large infrared excess and a flat UV slope (upper 
and lower panel of Fig. 7). Fainter bursts, which 
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Fig. 8. — Same as in Fig. 7 but for a random sam- 
ple of model bursts (crosses) , as described in the 
text. The filled pentagons are the same observa- 
tional data from Heckman et al. (1998). 



have lower gas column densities and are on average 
also more metal poor, quickly exhaust their gas 
content and evolve toward a low infrared excess 
and a negative UV slope, along the locus defined 
by observations (lower panel of Fig. 7). Conversely, 
brighter bursts, having larger gas column densities 
and higher metallicities, remain highly enshrouded 
by dust until, after a few tens of Myr, the dom- 
inant stellar population has become intrinsically 
old. Their UV continuum slopes always remain 
flat, at the beginning because of reddening and at 
later times because of age. 

In summary, the model bursts lie close to the 
region occupied by observed bursts in Fig. 7 as 
long as the stellar population is young, in the sense 
of the UV light being dominated by very massive 
stars. The position of bursts along the observa- 
tional locus is then determined mostly by the net 
dust opacity in the UV, in agreement with the 
interpretation of Meurer et al. (1995) and Heck- 
man et al. (1998). This in turn depends both on 
the initial gas mass, radius and metallicity of the 
burst, and on its evolutionary stage. 
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A detailed comparison with the Heckman et al. 
observations would require us to construct a mock 
catalogue of model starbursts obeying the same 
selection criteria as the observed sample. Unfor- 
tunately, the observational selection criteria are 
rather ill-defined. In addition, one of the selection 
criteria is the presence of strong HII region emis- 
sion lines, and the GRASIL code at present does 
not calculate these emission line properties. In- 
stead, we simply select starbursts with ages since 
the start of the burst less than t max , to account 
roughly for the effect that as soon as most of the 
massive stars have evolved away, the galaxy will 
no longer produce strong emission lines, and so no 
longer be classified as a burst in the observational 
sample. Fig. 8 shows the resulting distribution of 
points for the choice t max — 50Myr. The model 
starbursts are seen to follow similar relations to 
the observational sample. The results do not de- 
pend sensitively on the choice of t max . 

Several parameters may in principle affect the 
the spectral properties of a model starburst galaxy 
and therefore the location of our models in the 
above plots, but the most critical are the ratio 
between the star formation timescale and the dy- 
namical time (e~^ urst = T* huYSt / r hulgc ) , and the es- 
cape time (t csc ) for newly born stars to escape 
from their parental clouds. The former affects 
the bolometric luminosity, which is almost di- 
rectly proportional to the star formation rate, and 
thus inversely proportional to the star formation 
timescale. The latter affects the fraction of light 
absorbed inside clouds, and so may affect both the 
slope of the UV spectrum and the ratio between 
the IR and UV luminosities. The distribution of 
model points in Fig. 8 can therefore be used to 
constrain the values of e*burst and t csc . However, 
we found that changes in either of these parame- 
ters by a factor ~ 2 either way would only slightly 
worsen the match with observations. 

6.3. Infrared colours 

We now consider the infra-red and sub-mm 
colours of starbursts and normal galaxies. Fig- 
ure 9 shows the dependence of the mean IRAS 
colours on infra-red luminosity. This plot includes 
all model galaxies, both normal and starbursts. 
Their IRAS band luminosities are calculated by 
convolving the SEDs with the IRAS response func- 
tions. In calculating the mean colours, the models 
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Fig. 9. — Infra-red colours versus total infra-red 
luminosity. The three panels show the 12/25, 
12/60 and 60/100 ^m colours respectively, defined 
as the log of the ratio of vL v in each band. In 
each panel, the solid line with error bars shows the 
mean and dispersion from the models, while the 
dashed line with error bars shows the mean and er- 
ror on the mean measured by Soifer & Neugebauer 
(1991) from the IRAS Bright Galaxy Survey. 

are weighted by their number density and by a fac- 
tor L^y 2 (60^m) to account approximately for the 
volume within which a galaxy would be visible in 
a 60/im flux-limited sample. is the standard 
estimate of the total 8 — 1000/im IR luminosity 
from the luminosities in the four IRAS bandpasses 
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(Sanders & Mirabel 1996): 

L IR = 0.971/^(12) + 0.771/^(25) 

+ 0.931/1^(60) + 0.60^(100). (16) 

The model predictions are compared with the ob- 
served mean colours calculated by Soifer & Neuge- 
bauer (1991) from the IRAS bright galaxy sample 
(IRAS BGS, Soifer et al. 1989), which is a com- 
plete sample flux-limited at 60/xm. Models and 
observations, in particular the 12/60 /zm colour, 
show in general similar trends. 
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Fig. 10. — Average SED of model starbursts with 
Lir > 10 n /i~ 2 L Q (empty circles) compared to 
the observational data (crosses) for 14 luminous 
starbursts from Lisenfeld et al. (2000). (We ex- 
clude galaxies with only upper limits at 850^m.) 
All SEDs are normalized to the 60^m luminosity. 
The error bars show the dispersion in the models. 

Figure 10 shows the quite good agreement be- 
tween the average IR and sub-mm SED of model 
starburst galaxies, and the observed SEDs of lumi- 
nous infra-red galaxies from the sample of Lisen- 
feld et al. (2000). The limit L IR > lO n /i- 2 £ for 
the models has been chosen to approximately re- 
produce the selection for the Lisenfeld et al. sam- 
ple. Note that the dust opacity in our models de- 
creases as \- b with b « 2 for 100 < A < 1000/xm, 
while Lisenfeld et al. , by fitting optically thin sin- 
gle temperature models to the data at A > 60^iin, 



derived b values in the range 1.5-2. Our models 
demonstrate that the shallower slopes can instead 
be explained by the distribution of dust tempera- 
tures within each galaxy. 

6.4. Extinction in starburst galaxies 

An important problem in the study of star- 
forming galaxies is to determine the amount of 
attenuation of starlight by dust, especially in the 
UV. This bears directly on the determination of 
star formation rates in galaxies from their UV lu- 
minosities. For our own and a few nearby galax- 
ies, the extinction law of the dust can be mea- 
sured directly from observations of background 
stars, where the dust acts as a foreground screen. 
The differences found between the shapes of the 
extinction curves of the Galaxy, the Large Magel- 
lanic Cloud and the Small Magellanic Cloud be- 
low A <2600A (e.g. Fitzpatrick 1989) are often 
ascribed to the different metallicities in these sys- 
tems, covering the range Z ~ 0.1 — 1 Z@. Re- 
cently, Calzetti et al. (1994) (see also Calzctti 
1997, 1999) have analyzed the dust extinction in 
starburst galaxies. In this case, the measurement 
of the extinction is more complicated, since one 
measures the integrated light of the whole system, 
where stars and dust are mixed in a complex way. 
From the optical and UV spectra of a sample of 
UV-bright starbursts, Calzetti et al. derive an av- 
erage attenuation law characterized by a shallower 
far-UV slope than that of the Milky Way extinc- 
tion law, and by the absence of the 2175 A feature. 
This is at first sight quite surprising, because the 
metallicities of these galaxies are mostly similar to 
that of the Milky Way, and so they might be ex- 
pected to have similar dust properties. The ques- 
tion is then to what degree the differences between 
the starburst attenuation law and the Milky Way 
extinction law are due to the geometry of the stars 
and dust, and to what degree they can only be ex- 
plained by differences in dust properties. 

Figure 11 compares the average attenuation 
curves for galaxies from our model with the em- 
pirical "attenuation law" obtained for starbursts 
by Calzetti et al. (1999). The attenuation A\ for 
the models is defined as the difference in magni- 
tudes of the stellar luminosity L\ of a galaxy with 
and without dust, and is normalized to the colour 
excess E(B — V) = As — Ay of the stars to give 
an attenuation "law" fc(A) = A\/E(B — V), equiv- 
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Fig. 11. — The average dust attenuation curves 
for starlight in different classes of galaxies (normal 
and starburst, SB) in the model compared with 
the average Milky Way extinction law (solid line) 
and with the Calzetti "attenuation law" (filled cir- 
cles, Calzetti ct al. (1999), with R' v = 4.05). The 
attenuation curves are all normalized to the net 
reddening of the stellar population, E(B — V). 
The error bars show the dispersion of the mod- 
els around the mean attenuation curve. 

alent to the definition of Calzetti et al. . As de- 
scribed in §3.2, the dust properties we adopt imply 
an extinction law characterized by a distinct 2175 
A feature produced by graphite grains, and well 
matching the average Milky Way extinction curve. 
The model extinction law (solid line in Fig. 11) is 
the attenuation law that would be measured if all 
the dust were in a foreground screen in front of the 
stars and no scattered light reached the observer. 
This geometry is clearly not realistic as applied to 
the integrated light from galaxies. In our models, 
we have instead a complex and wavelength depen- 
dent geometry, where the UV emitting stars are 
heavily embedded inside molecular clouds, while 
the older stars, mainly emitting in the optical and 
near infrared, are well mixed with the diffuse in- 
terstellar medium. 

Figure 11 shows average attenuation curves for 
3 classes of model galaxies: (a) normal galax- 
ies with E(B - V) > 0.05; (b) starbursts with 



5 x 10 8 < L IR < 5 x 10 w h- 2 L Q ; and (c) star- 
bursts with L IR > 5 x 10 10 /i~ 2 L Q . The star- 
burst models are all chosen to have ages < 50Myr 
since the start of the burst, as discussed in §6.2. 
Sample (b) corresponds roughly to the galaxies for 
which Calzetti et al. measured their attenuation 
law. The model attenuation law depends signifi- 
cantly on the sample, but all 3 classes show a weak 
or completely absent 2175A feature. In particular, 
the predicted attenuation curve for the lower lu- 
minosity starbursts is remarkably close to the em- 
pirical "Calzetti law". This result is an entirely 
geometrical effect, and did not require us to as- 
sume for starbursts dust properties different from 
those of the Galaxy. This conclusion is contrary to 
that of Gordon ct al. (1997), who argued that the 
observed shape is only produced with dust that 
lacks the 2 175 A feature in its extinction curve. 
The reason is presumably that Gordon et al. only 
considered clumping of dust, not of stars, and as- 
sumed a spatial distribution for stars independent 
of stellar age. Our results follow naturally from 
the assumption that stars are born inside dense 
dust clouds and gradually escape. 




4000 8000 
X[A] 



Fig. 12. — Predicted attenuation curves for two 
normal galaxies (left panels) and two starbursts 
(right panels) from the models. 

To further illustrate the importance of geomet- 
rical effects in determining the attenuation law, 
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we show in Fig. 12 the attenuation laws of 2 nor- 
mal and 2 starburst model galaxies. The global 
attenuation law, and the separate contributions 
from the molecular clouds and diffuse dust, are 
shown in each case, normalized to the colour ex- 
cess E(B — V) produced by that dust component. 
The global (g), molecular cloud (MC) and diffuse 
dust (d) contributions are related by 



(A x /E) MC Emc + (A x /E) d E d 
E„ 



(17) 



In the far-UV, including the spectral region 
around the 2 175 A feature, the global attenua- 
tion in the models is strongly contributed, or even 
dominated, by the MCs. The shape of the attenu- 
ation curve there has little to do with the optical 
properties of grains, because our MCs usually have 
such large optical depths that the UV light from 
stars inside the clouds is completely absorbed. 
The wavelength dependence of the attenuation 
law of the MC component instead arises from the 
fact that the fraction of the light produced by very 
young stars increases with decreasing wavelength, 
and at the same time, the fraction of stars which 
are inside clouds increases with decreasing age, 
as given by eqn. (11). The additional attenuation 
arising in the cirrus component can sometimes im- 
print a weak 2175A feature, but this is not the case 
for the starbursts, where the primary UV stellar 
light is dominated by very young populations. 

7. Galaxy Luminosity Function 
7.1. Method 

The luminosity function of galaxies at different 
wavelengths is a basic property of the galaxy pop- 
ulation which a galaxy formation model should ex- 
plain. We calculate the galaxy luminosity function 
at different wavelengths by combining the model 
SEDs with the weights for the individual galax- 
ies (Section 4). For the normal galaxy sample we 
have, for the number density of galaxies per lni 
at some wavelength A 
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(18) 



where rii is the number density for the ith galaxy, 
Li is its luminosity at wavelength A, the centre 
of the bin is at L and its width is A(lnL). For 



the burst galaxy sample, we have to sum over the 
burst phase j also, giving 
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| lnL„-lnL|<iA(lnL) 



where n^ is the number density for galaxy i at 
evolutionary phase j, and its the luminosity 
at that phase. 

Galaxy luminosity functions are measured in 
specific bands defined by a filter+instrument re- 
sponse function, e.g. the standard B or K bands, 
or the IRAS bands. Thus we convolve the model 
SEDs with the appropriate response function 
to calculate the luminosity L v in that band. 
We use absolute magnitudes on the AB system, 
M AB = -2.51og 10 (L„/4.345 x KPergs^Hz" 1 ). 
The model luminosity functions have statistical 
uncertainties due to the finite size of the model 
galaxy catalogue. We estimate these statistical 
errors by bootstrap resampling of the catalogue. 

7.2. Optical and Near Infra-Red 

At optical and near-IR wavelengths the emis- 
sion is mostly from older stars and the effects of 
dust obscuration are generally modest. Figure 13 
shows the local luminosity function in the B-band 
(A = 0.44^m), compared to the observed lumi- 
nosity function measured from the ESP redshift 
survey by Zucca et al. (1997). The predicted lu- 
minosity function agrees well with the observed 
one, except at the highest luminosities. Extinction 
by dust makes galaxies around 0.6mag fainter on 
average, for bright (L > L+) galaxies. Galaxies 
which have had recent bursts (i.e. in the last 1/20 
of the age of the universe, 0.7 Gyr) do not dom- 
inate the luminosity function at any luminosity, 
when the effects of dust are included. 

As described in Cole et al. (2000), the B- 
band luminosity function is used as one of the 
primary observational constraints for setting the 
parameters in the GALFORM model, in particu- 
lar, the parameters ahot an d Vhot controlling feed- 
back, and the parameter T which sets the fraction 
of brown dwarfs in the IMF. The good agreement 
with the observed B-band luminosity function is 
therefore not a surprise, but it was not guaran- 
teed, since the stellar population and dust mod- 
els used in Cole et al. (2000) are not identical 
to those used here. Cole et al. used the stellar 
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Fig. 13. — Galaxy luminosity functions in the optical and near-IR. The left panel shows the luminosity 
function in the B-band (0.44/im), compared to observational data from Zucca et al. (1997). The right 
panel shows the luminosity function in the K-band (2.2/im), compared to observational data from Gardner 
et al. (1997). In both panels, the solid line shows the total luminosity function including the effects of 
dust, while the long-dashed and short-dashed lines show the contributions to this from galaxies with and 
without recent bursts, and the dotted line shows the luminosity function without dust. The error bars on the 
model luminosity functions are the statistical errors resulting from the finite sample size, calculated using 
bootstrap resampling. The observational data have been converted to AB magnitudes assuming B = Bj +0.2 
(to convert from measured to standard magnitudes) and Bab = B — 0.12, and Kab = K + 1.87 respectively. 



population models of Bruzual & Chariot (2000), 
and calculated the effects of dust using the mod- 
els of Ferrara et al. (1999). The stellar popula- 
tion model in GRASIL is based on similar stellar 
evolution tracks and spectra, but the treatment 
of dust extinction is significantly different. The 
Ferrara et al. models assume that stars and dust 
are smoothly distributed, while in GRASIL a frac- 
tion of the dust is in clouds, and young stars are 
confined to these clouds. The B-band luminos- 
ity functions, both with and without dust, calcu- 
lated by GALFORM using the Bruzual & Char- 
lot (2000) and Ferrara et al. (1999) models agree 
very well with those computed using the GRASIL 
stellar population+dust model, demonstrating the 
consistency of the procedure of using the galaxy 
formation parameters derived in Cole et al. in com- 
bination with the GRASIL model. The effects of 
dust computed using the two models are quite sim- 
ilar in the B-band, in spite of the differences in the 
star and dust geometry. This is because most of 
the B-band light is produced by stars which are old 
enough to have escaped from the clouds in which 



they formed, so in GRASIL the attenuation is due 
mostly to the diffuse component of the dust, which 
is modelled in a similar way to that in the Ferrara 
et al. models. 

Figure 13 also shows the model and observed lu- 
minosity functions in the K-band. In this case, the 
effects of dust are very small, so the comparison 
is essentially independent of assumptions about 
dust. Again, the model agrees well with obser- 
vations over most of the luminosity range, as was 
also found by Cole et al. (2000). The contribution 
of galaxies with recent bursts is very small at all 
luminosities. 

7.3. Far Ultra- Violet 

In Figure 14 we compare the predicted lumi- 
nosity function in the far-UV (A = 0.2/xm) with 
that measured by Sullivan et al. (2000) from a 
UV-selected redshift survey, based on FOCA in- 
strument fluxes. This comparison has not previ- 
ously been made for any semi-analytical galaxy 
formation models. The effect of dust are much 
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Fig. 14. — Far-UV (0.2/rni) luminosity function, (a) The left panel shows the results for our standard 
parameters. The solid line shows the model prediction for the LF at z — including dust, while the 
long-dashed and short-dashed lines show the contributions to this from galaxies with and without recent 
bursts. The dotted line shows the model luminosity function obtained if absorption by dust is neglected. 
The symbols with error bars show the observed 0.2/xm luminosity function (without any dust correction) 
measured from a UV-selected redshift survey by Sullivan et al. (2000). FOCA magnitudes are converted 
to the AB system, assuming itiab — ^2000 + 2.26. (b) The right panel shows the effects of varying t esc 
and rburst- The solid line is for our standard model (t csc — 2 and lOMyr in normal galaxies and bursts 
respectively, and rburst Abuigc = 0.1), including dust. The short-dashed line shows the effect of increasing 
t esc (burst) to 30Myr, and the short-dash-long-dash line the effect of increasing t esc (normal) to 5Myr. The 
long-dashed line shows the effect of increasing rb U rst/f bulge to 0.5. The dotted line shows the LF without 
dust, which is the same in each of these cases. The observational data are as in the left panel. 



larger than in the optical, as one would expect. In 
this case, the effects of the more realistic geometry 
for the stars and dust assumed by GRASIL com- 
pared to the Ferrara et al. models (clumpy rather 
than smooth distributions for the stars and dust) 
are significant. The stars that produce most of 
the UV light spend a large fraction of their life- 
times in the molecular clouds where they form, so 
the mean extinction is larger than in the case of 
a smoothly distributed dust component with the 
same total dust mass. Bursting and non-bursting 
galaxies contribute roughly equally at the highest 
luminosities. This result is however sensitive to 
the details of how bursts are modelled, since this 
determines what small fraction of the UV light es- 
capes from currently or recently bursting galaxies. 
When we compare our model LF including extinc- 
tion with the directly observed LF, uncorrected for 
extinction, we find reasonable agreement at lower 
luminosities, but at high luminosities, the model 



LF is somewhat lower than the observed one. This 
might be partly an effect of evolution in the obser- 
vational sample, which covers a significant redshift 
range (z < 0.5), but it might also be that the UV 
extinction is over-estimated in the model. 

Figure 14 shows also the effect of changing the 
burst radius rburst and the timescale t esc for stars 
to escape from clouds. Increasing rburst /''bulge 
from 0.1 to 0.5 reduces the optical depth in the 
diffuse component during bursts, allowing more of 
the UV light from bursts to escape, and increas- 
ing the LF at the highest luminosities. Increasing 
t osc in bursts from lOMyr to 30Myr has negligi- 
ble effect on the total UV LF. Increasing t csc in 
normal galaxies from 2Myr to 5Myr slightly low- 
ers the amplitude of the luminosity function at the 
bright end. 
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Fig. 15. — The 60^m far-IR luminosity function, (a) The left panel shows the results for our standard 
parameters. The solid line shows the total model luminosity function, while the long-dashed and short- 
dashed lines show the contributions to this from galaxies with and without recent bursts. The points with 
error bars show the observed 60^m luminosity functions measured from galaxies observed by IRAS, by 
Saunders et al. (1990) (filled circles) and Soifer & Neugebauer (1991) (open circles), (b) The right panel 
shows the effects of varying the burst Or eSCcip6 tim6SCcil6SI ^*burst = T*burst/Tbuigc = 1 (dotted), 2 (solid), 4 
(dashed). The long dashed curve has e* burst = ^ ( our standard value) but t csc = 2 Myr instead of 5 Myr. 



7.4. Mid and Far Infra-Red 

In the mid- and far-infrared, the luminosity of 
galaxies is dominated by re-emission from dust. 
Using the GRASIL code, we can now directly pre- 
dict the far-IR luminosities of galaxies from our 
galaxy formation model, and compare with obser- 
vations. The luminosity functions of galaxies at 
12, 25, 60 and lOO^m have been measured using 
IRAS data. The best determination is at 60^im, 
where IRAS was most sensitive. Figure 15 shows 
that the predicted luminosity function agrees ex- 
tremely well with that observed by Saunders et al. 
(1990) and Soifer & Neugebauer (1991), except at 
very low luminosities, where the measured LF is 
fairly uncertain. Above ^£^(60) > 3 x 10 10 /i~ 2 £ Q , 
the model LF is dominated by galaxies undergo- 
ing bursts triggered by mergers. This is in agree- 
ment with observations of ultra-luminous IRAS 
galaxies, which are all identified as recent mergers 
based on their optical morphology (e.g. Sanders 
& Mirabel 1996). 

The right panel of Figure 15 shows the ef- 
fect on the 60/im LF of varying the parameter 
e*burst, which relates the star formation timescale 



in bursts to the dynamical time of the bulge (equa- 
tion 6). Unlike the other parameters in the GAL- 
FORM model, Cole et al. (2000) did not try 
to choose a best-fit value, because the observa- 
tional data in the optical and near-IR that they 
compared with were not sensitive to its value. 
(The Cole et al. results were calculated assuming 
T*burst = 0.) However, the far-IR LF is sensitive 
to this and thus constrains the burst timescale for 
the most luminous galaxies. Figure 15 shows pre- 
dictions for e*burst = 1, 0.5, 0.25, corresponding 
to r^burst/Tbuigo = 1,2,4 respectively. Increasing 
e*burst means bursts are more luminous, but last 
for a shorter time, and so have a lower number 
density. This trend is seen at the high-luminosity 
end of the 60/zm LF, which is dominated by burst- 
ing galaxies. A value e*burst = 2 seems to fit some- 
what better than higher or lower values, so we 
adopt this as our standard value. Also shown in 
the same panel is the somewhat better fit obtained 
setting t csc — 5 Myr. However, as explained in 
§ 3.3, our adopted standard value 2 Myr is favored 
by stellar evolution timescale argument and by the 
UV LF. Increasing t osc in bursts from lOMyr to 
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30Myr has negligible effect on the LF. 
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Fig. 16. — The luminosity functions in the 12, 
25 and lOO^im IRAS bands, compared to obser- 
vational data from Soifer & Neugebauer (1991). 
The line types and symbols are as in Fig. 15. 

The luminosity functions at 12, 25 and 100 /im, 
are compared with the observational data from 
Soifer & Neugebauer (1991) in Figure 16. The 
predicted luminosity function agrees well with the 
measured one in each case. 

8. Star formation rate indicators 

Here we examine the accuracy of several SFR 
indicators based on continuum UV or IR luminosi- 
ties (reviewed by e.g. Kennicutt 1999). 

The luminosity £„(2800) at 2800 A has been 
extensively used to estimate SFRs of high-redshift 
galaxies and to investigate the evolution of the cos- 
mic SFR density (e.g. Lilly et al. 1996; Conolly 
et al. 1997). In the top panel of Fig. 17, we plot 
the SFR against L„(2800) for the model galaxies, 
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Fig. 17. — Star formation rate versus UV and IR 
luminosities for galaxies from our model. Only 
models corresponding to spirals and luminous 
starburst galaxies are plotted. See text for details. 

including the effects of extinction. Only models 
corresponding to spiral galaxies (B/T < 0.5) and 
luminous starburst galaxies (L/ij > 10 10 h~ 2 L Q ) 
are shown. At higher luminosities and SFRs, 
£,,(2800) > 3 x lO^ergs^Hz -1 , the models 
in the absence of dust follow a linear relation 
between SFR and L„, with a rather small dis- 
persion, as would be expected if the UV lu- 
minosity is dominated by young stars and the 
recent SFR has been approximately constant. 
This linear relation, S'Fi?/(M Q yr~ 1 ) = 8.5 x 
10- 29 J L^(2800l)/(ergs- 1 Hz" 1 ), is indicated by 
the solid line, and its extrapolation to lower lumi- 
nosity is shown by the dotted line. Dust extinction 
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shifts points to the left of this line. At lower lu- 
minosities, the effects of dust extinction are very 
small, because of the low gas contents of the galax- 
ies. On the other hand, because the SFRs are so 
small, the 2800A light has a significant contri- 
bution from post-AGB stars and old metal poor 
populations, and this causes the locus of points 
to bend to the right of the linear relation. The 
galaxies with very high mctallicitics (Z > 0.1) 
have very low gas fractions. 

The figure also shows as a dashed line the linear 
relation between L I/ (2800) and SFR obtained by 
Kennicutt (1999), using stellar population mod- 
els for a Salpeter IMF, and assuming a constant 
SFR for the last 10 8 yr. The SFR/L V ratio in 
our models without dust (solid line) is about 40% 
lower than Kennicutt's value, but this difference is 
entirely due to the different IMF we adopt (equa- 
tion (9)). 

Perhaps the most striking feature of this plot 
is that the starburst models are offset by more 
than an order of magnitude from the average rela- 
tion holding for the normal spirals, because of the 
large UV extinctions in the starbursts. Further- 
more, their dispersion in SFR at a given luminos- 
ity is also quite large. Thus, the 2800A luminosity 
with no dust correction performs rather poorly as 
a quantitative SFR indicator, both for very high 
SFRs (because of extinction) and for very low ones 
(because of the light from older stars) . 

The middle panel of Fig. 17 depicts the relation 
between the star formation rate and (eq. 16), 
the estimated 8 — 1000/im luminosity based on the 
IRAS fluxes. The solid line represents the rela- 
tion derived by Kennicutt (1998) for starbursts, 
SFR/iMoyr' 1 ) = 4.5 x 10" 44 W^gs" 1 ), by 
assuming that the bolometric output in a continu- 
ous burst of age between 10-100 Myr is completely 
reprocessed by dust, again for a Salpeter IMF. The 
Kennicutt relation is seen to fit our model galax- 
ies quite well at all luminosities (normal spirals 
as well as starbursts), even though we assume a 
different IMF from Kennicutt. 

The luminosity in the ISO 15^m band has also 
been proposed as an approximate SFR indica- 
tor. The bottom panel of Fig. 17 shows the SFR 
vs -L„(15/im) for our galaxies. The line plot- 
ted is our best linear fit to the model points, 
SFR/iMQyi- 1 ) = 5.6xl0- 3 ° J L I ,(15^m)/(ergs- 1 Hz" 
Despite the important contribution from PAH 



bands, the correlation of the SFR with L„(15/im) 
in our models is still fairly good. 

9. Summary and Conclusions 

We have combined an ab initio model of galaxy 
formation (GALFORM, §2, Cole et al. 2000), with 
an ab initio model for stellar emission and dust 
emission and absorption in galaxies (GRASIL, §3, 
Silva et al. 1998). Both models are state-of-the 
art. We are able to predict, in the context of the 
cold dark matter cosmology, the luminosities and 
spectral energy distributions from the UV to the 
sub- mm for the whole galaxy population, and how 
these change with cosmic epoch. Here we focused 
on a wide range of spectrophotometric properties 
of present-day galaxies, from the UV to the sub- 
mm, for a flat low-density cosmology (fio = 0.3, 
A = 0.7) with a CDM spectrum of density fluctu- 
ations. The model is remarkably successful in ex- 
plaining the UV, optical and IR spectrophotomet- 
ric and extinction properties of galaxies in the lo- 
cal universe. Future papers will investigate galaxy 
evolution in the UV, optical, IR and sub-mm out 
to high redshift. 

Dust plays a dominant role in starburst galax- 
ies, where star formation proceeds in the central 
regions of a galaxy on a short timcscale. We did 
not previously make any detailed comparison of 
the properties of the starbursts predicted by semi- 
analytical galaxy formation models with observa- 
tional data. Here we have shown that these prop- 
erties are nicely reproduced. 

Our model predicts an average dust attenuation 
law for starburst galaxies that agrees remarkably 
well with the empirical law found by Calzetti et 
al. (1999), although with a significant dispersion 
around the mean. In particular, the 2175A bump 
is absent when the net attenuation of the galaxy 
light is considered. This is entirely an effect of the 
geometry of stars and dust in our model, and has 
nothing to do with the optical properties of dust 
grains. Indeed, our dust mixture would reproduce 
the average Milky Way extinction curve (with a 
strong 2 175 A feature), if it were arranged in a 
foreground screen geometry. The absence of this 
feature in the attenuation curves of model star- 
bursts is because in that case the dust attenua- 
tion is dominated by molecular clouds, with the 
shape controlled by the gradual escape of young 
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stars from the clouds. 

The starburst galaxies are predicted to dom- 
inate the bright end of the luminosity function 
in terms of bolometric luminosity, but because of 
the large extinctions in these objects, they do not 
make a dominant contribution to the bright end 
of the luminosity function in cither the UV, op- 
tical or near-IR once dust effects are included. 
At these wavelengths, the luminosity function is 
dominated by normal spiral and elliptical galax- 
ies. However, the starbursts completely dominate 
the bright end of the luminosity function in the 
mid- and far-IR (10-100/im), at total luminosi- 
ties Ljr > 10 n /i~ 2 L Q . Overall, the luminosity 
function predictions from the far-UV to the far- 
IR are a remarkable success for the model, since 
the dust contents and galaxy radii are predicted a 
priori, and the only significant adjustable param- 
eter in the comparison was the ratio of the burst 
timescale to the bulge dynamical time, which was 
chosen to fit the bright end of the luminosity func- 
tion at 60yum. 

As expected, our models show that the UV con- 
tinuum is in general a poor star formation indica- 
tor, both because of the large variations in the 
amount of extinction, and also because of the con- 
tribution from old stellar populations in the mid- 
UV (~ 3000A) in more quiescent galaxies. The 
infrared luminosity is a much more reliable SFR 
indicator. 

The parameters values adopted here for the 
GALFORM model are those chosen previously by 
Cole et al. (2000) to fit the properties of the local 
galaxy population in the optical and near-infrared, 
apart the new ones specifying the timescales and 
radii of bursts. Cole et al. also discuss the effects 
of varying the 'old' GALFORM paramaters. The 
purpose of this paper was to present the effects of 
including dust in a fixed galaxy formation model. 
The treatment of dust reprocessing with GRASIL 
requires some additional parameters to be set, but 
opens the possibility to test semi-analytical mod- 
els against the wealth of IR and sub-mm obser- 
vations already available or planned for the near 
future. The adopted values for these parameters 
have also been guided by the results of Silva et 
al. (1998), who used GRASIL to reproduce de- 
tailed SEDs of several local normal and starburst 
galaxies. Some of the effects of variations of these 
parameters on SEDs, LFs and starburst properties 



have been discussed in the relevant sections of this 
paper. A more systematic parameter study will 
be included in future investigations. Among the 
newly introduced parameters, probably the most 
important ones are e+burst and t CSCl describing the 
timescales of bursts and the time for young stars to 
escape from their parent molecular clouds. They 
have significant effects on (and are constrained by) 
the IR and UV luminosity functions respectively. 

In conclusion, this paper is a stepping stone 
for future work, which will apply the same mod- 
els to galaxies at high redshift. Now that semi- 
analytic models can be effectively compared with 
infrared and sub- millimetre observations, as well 
as UV and optical data, they can be used to work 
towards an observationally and theoretically con- 
sistent picture for the history of galaxy formation 
and star formation in the universe. 
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